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Abstract 

A mesoscopic continuum model is employed to analyse the transport mechanisms and structure 
formation during the redistribution stage of deposition experiments where organic molecules are 
deposited on a solid substrate with periodic stripe-like wettability patterns. Transversally invariant 
ridges located on the more wettable stripes are identified as very important transient states and 
their linear stability is analysed. It is found that there exist two different instability modes that 
result (i) at large ridge volume in the formation of bulges that spill from the more wettable stripes 
onto the less wettable bare substrate and (ii) at small ridge volume in the formation of small 
droplets located on the more wettable stripes. These predictions are confirmed by direct numerical 
simulations of the fully nonlinear evolution equation for two-dimensional substrates. In addition, 
the influence of different transport mechanisms during redistribution is investigated focusing on 
the cases of convective transport with no-slip at the substrate, transport via diffusion in the film 
bulk and via diffusion at the film surface. In particular, it is shown that the transport process does 
neither influence the linear stability thresholds nor the sequence of morphologies observed in the 
time simulation, but only the ratio of the time scales of the different process phases. 
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I. INTRODUCTION 


Many coating and surface growth processes through which a solid substrate is covered 
by homogeneous or structured layers of the same or other materials combine deposition 
and re-distribution stages that might occur successively or in parallel. Prominent examples 
are various homo- and heteroepitaxial surface growth processes where the re-distribution 
of material on the substrate occurs through diffusion of deposited atoms or molecules on 
the surface of the substrate [ETB06, EDM13]. Often, the material is first deposited, e.g., 
via vapor deposition [WZJN01], molecular beam deposition [Sie97] or (pulsed) laser depo¬ 
sition [MEvdB05]. Then the adsorbed molecules diffuse along the surface and form two- 
dimensional aggregates, complete mono- or multilayers, and three-dimensional nanoscale 
structures like pyramids, holes and mounds - the latter are sometimes called quantum dots 
and wells [GDV03, KE10]. 

Another example is spin-coating where first a drop of liquid is deposited on the sub¬ 
strate, before convective motion caused by fast spinning of the substrate spreads the liquid 
into a thin film [WHDOO, SR04], In the case of solutions or suspensions the spin coating 
process is often accompanied by solvent evaporation and the final coating is a homogeneous 
or structured layer of the dried-in solute [RBD91, TMP98, MPW11]. Another example are 
dip-coating and related processes where a liquid is transferred from a bath or other reser¬ 
voir onto the substrate [MOY03, Thil4b]. These examples all involve a wet stage where 
hydrodynamic flows caused by external forces, wettability and capillarity are important for 
the redistribution of the material [WR04]. The flows during the wet stage may be described 
employing, e.g., asymptotic models as obtained via long-wave expansions from the basic 
equations of hydrodynamics [ODB97]. 

Many of the studied systems involve homogeneous substrates. However, there is also a 
growing number of experimental and theoretical studies that investigate the use of hetero¬ 
geneous substrates to control the structure of deposits. Examples include recent deposition 
experiments with organic molecules performed by Wang et al. [WDW + 11, WC12] on silicon 
oxide substrates with gold stripes, the study of dip-coating for chemically micropatterned 
surfaces [DTD + 00] or dewetting of thin silicon films on nano-patterns formed by electron 
beam lithography [BAR + 13]. For the growth of quantum dots on heterogeneous substrates 
see Ref. LWY14. 
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Static morphologies emerging from capillary and wettability influences on heterogeneous 
substrates of various geometries are already studied in depth, however, studies of the in¬ 
volved dynamics are less frequent: Morphological transitions of transversally invariant liq¬ 
uid ridges on two-dimensional striped substrates (i.e., mathematically equivalent to drops 
on one-dimensional heterogeneous substrates) were studied via minimization of macroscopic 
interfacial free energies in Ref. LL98 and through bifurcation studies employing mesoscopic 
free energies (interface hamiltonians) and time simulations of gradient dynamics on these en¬ 
ergies (equivalent to long-wave hydrodynamic models) in Refs. BKTB02, TBBB03. In these 
studies nonvolatile liquids are considered, i.e., liquid volume and heterogeneity properties 
are the main control parameters. Related steady (one-dimensional) results are obtained in 
studies of morphological transitions of wetting films on striped substrates [BDP99, BDOO]. 
There, the main control parameters are temperature and chemical potential, i.e., the volume 
of deposited liquid is a dependent quantity in this grand canonical studies. 

Fully two-dimensional situations are also studied: Ref. TBBB03 investigates the linear 
instability modes and their time scales for transversally invariant liquid ridges on homo¬ 
geneous and heterogeneous (striped) substrates. In particular, they study transversal in¬ 
stabilities (Plateau-Raylcigh instability) and their coupling to coarsening instabilities while 
Ref. BKHT11 also considers the coupling to depinning modes under lateral driving and the 
bifurcation structure for steady two-dimensional states (i.e., height profiles depend on both 
substrate dimensions). The Plateau-Rayleigh instability is also considered in Ref. MRD08 to¬ 
gether with the effect of a body force along the ridge. Sharma et al. present time simulations 
of dewetting dynamics for one-dimensional heterogeneous substrates and two-dimensional 
substrates with a less wettable square patch [KKSOO], with stripes [KS01] and other two- 
dimensional wettability patterns [KS03]. Droplet spreading on patterned substrates is also 
considered [VSK11], 

Most of the mentioned studies that involve two-dimensional substrates employ meso¬ 
scopic long-wave models (i.e., small gradient expansions) while Refs. BL02, BKL05 deter¬ 
mine steady surface profiles and their stability on the basis of macroscopic interfacial free 
energies (note, that this does not allow for a calculation of the time scales of the instabilities 
as they result from a balance of dissipation mechanisms and decrease in free energy). Similar 
results for morphological transitions on substrates with micro-grooves are compared with 
experiments in Ref. SBK + 05. Other experiments concern the pattern-directed dewetting of 
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ultrathin polymer films that results, e.g., in rows of drops [SFD + 02], Macroscopic exper¬ 
iments with water on chemically heterogeneous substrates are reported in Ref. GHLL99. 
Part of these works are reviewed in Ref. HBS08. 

The present study is directly motivated by deposition experiments with organic mol- 
cules performed by Wang et al. [WDW + 11, WC12], Several types of light-emitting or¬ 
ganic molecules (e.g., diferrocene (DiFc); 1,6- Bis(2-hydroxyphenol)pyridinel boron bis (4- 
n-butyl-phenyl)-phenyleneamine ((dppy)BTPA); N,N’-Di[(N-(3,6-di-tert-butyl-carbazyl))- 
n-decyl] qninacridone (DtCDQA); N,N’-bis-(l-naphyl)-N,N’-diphenyl-l,l , -biphenyl-4,4’- 
diamine (NPB)) are deposited through vapour deposition onto Si0 2 substrates that are 
prestructured with parallel gold stripes of various width. In such vapour deposition exper¬ 
iments the ongoing processes often resemble epitaxial growth of inorganic materials: the 
molecules diffuse on the substrate, nucleate at favourable sites and form growing islands. 
The prepatterning allows for a control of the nucleation and growth [WC12], 

Although the deposition and redistribution stages often occur successively, deposition may 
also still continue when the material is already being redistributed via diffusion or convective 
motion. Note that in particular for highly mobile deposit material, the distinction between 
the two transport processes, diffusion and convection, is not sharp when one deals with 
deposit amounts that allow at least locally for multilayer structures. Then one is already 
able to define a velocity profile within the film . In consequence, depending on the particular 
material properties and annealing parameters one may expect that depending on the location 
on the substrate transport is dominated by diffusion in a surface layer of the deposit (as for 
solid films discussed, e.g., in [SS86, MPLS10]), by diffusion of the bulk layer (as in dynamical 
density functional theory for diffusion of an adsorbate layer [ART 10]), by convection with 
strong slip at the substrate (as for polymeric layers [MWW05]), or by convection without 
slip (as for the majority of liquid layers [ODB97, ThilO]). In the long-wave models for the 
time evolution for surface profiles of the deposition layers or, more generally, the adsorption 
at the substrate these different transport mechanisms are related to power-law mobilities of 
different powers (zero, one, two and three, respectively). A central aim of the present work is 
to clarify in which way a distinction of redistribution via the different transport mechanisms 
is important for the observed phenomena. 

This is an important question as, e.g., in Refs. WDW + 11, WC12 it is discussed that one 
of the employed molecules, DtCDQA, shows a liquid-like behavior while the other molecules 
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(as e.g. NPB) behaves like a solid. In the latter, solid case, the height of the rather flat 
deposits always decreases with their width. In contrast, in the former (liquid ) case, with 
increasing amounts of deposited molecules the DtCDQA assembles into droplets, elongated 
drops, molecule stripes, cylindrical ridges and ridges with bulges [WDW + 11]. The height 
of the ridges increases with their width, an effect attributed to capillarity (constant contact 
angles) [WC12], Examples of structures obtained in these experiments are reproduced in 
Figs. 1 and 2. 



FIG. 1: AFM image of DtCDQA molecules deposited on a SiC >2 substrate patterned with Au 
stripes. The width of the Au stripes is increasing from right to left from 0.3 pm to 2.3 pm. The 
mean film thickness is 50 nm. Since the diameter of one molecule is about 1 nm, this corresponds 
to roughly 50 monolayers. Reprinted with permission from WC12. Copyright 2015 American 
Chemical Society. 

Note, that the description as ’liquid-like behaviour’ in Ref. WDL + 13 is mainly based on 
the observed shape of the deposits (ridges with bulges) but not on an actual observation 
of the transport process or velocity profile across the deposit layer. Therefore, one might 
prefer to describe the behaviour as dominated by interfacial tensions, i.e., capillarity. Other 
capillarity-dominated morphologies are obtained for DtCDQA on substrates with other pat¬ 
terns, e.g., with a cross-bar structure. Aspects of these experiments have already been 
described via kinetic Monte Carlo simulations on a lattice [LMW + 12], in particular, the 
transitions between growth that is localized on the gold stripes, the development of bulges 
that partly cover the bare substrate and island formation everywhere on the susbtrate - that 
occur in dependence of the various interaction parameters. 
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FIG. 2: Different growth regimes during the deposition of DtCDQA molecules on Si02 substrates 
patterned with Au stripes of alternating width. Panels a) and b) show AFM images of 3nm 
and lOnrn DtCDQA, respectively. Panel c) shows sketches of the different growth regimes. For 
small amounts of deposited molecules, small droplets form on the Au lines (I). With increasing 
amount, larger droplets result until droplet edges reach the borders of the stripe (II) which become 
elongated (III). Finally, the whole stripe is covered with a cylindrical ridge (V). At a certain amount 
of molecules, bulges form that also cover part of the Si02 substrate. Picture reproduced from the 
Supplementary information section of Ref. WDW + 11. Copyright (c) 2015 WILEY-VCH Verlag 
GmbH Sz Co. KGaA. 

Transformations between liquid ridges with and without bulges are also observed in 
macroscopic experiments with water [GHLL99]. There, bulges form when the contact angle 
for contact lines pinned at the step in wettability exceeds 90°. As the equilibrium contact 
angles in the experiments with DtCDQA on Au-striped SiC >2 are below 22°, the macroscopic 
theory does not directly apply to the experimental findings. 

Here, we employ a mesoscopic continuum description to describe the deposition exper¬ 
iments. In particular, we are interested in the stability of transversally invariant ridges of 
material located on the more wettable gold stripes and, in general, the role of such ridge 
states in the course of the time evolution from a homogeneous deposited film of molecules 
towards a final bulge or drop geometry. We also consider the question of how the dominant 
redistribution mechanism (diffusion or convection) influences the evolution pathway. As 
our interest is in the static states and the dynamic behaviour and because contact angles 
are small, we employ a gradient dynamics in small-slope approximation (i.e., a thin him or 
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long-wave model) [ODB97, CM09, ThilO], which describes the temporal evolution of the 
film height profile of a thin layer of material as driven by wettability and capillarity. The 
heterogeneity of the substrate is modelled as a chemical heterogeneity that only affects the 
wettability similar to Refs. KKSOO, TBBB03, KS03, BKHT11 where a number of different 
wetting energies (and therefore Derjaguin pressures) are used. 

The focus of the analysis presented here lies on the different instabilities of a liquid ridge 
that result either in the formation of large bulges that spill from the more wettable stripes 
onto the less wettable bare substrate or in the formation of small droplets on the more 
wettable stripes. We perform a linear stability analysis of steady ridges using continuation 
techniques [KOGV07, DWC + 14] as outlined in TBBB03 and presented in tutorial form in 
Ref.Thil4a. In particular, we analyze the influence of the wettability contrast, the amount 
of deposited material (mean film thickness), and the geometry of the stripe pattern on the 
stability of liquid ridges. Particular attention is given to an analysis of the influence of the 
main transport mechanism in the redistribution stage. 

Our analysis allows us to identify and characterize the different growth regimes that 
are experimentally observed in Ref. WDW + 11 (here reproduced in Fig. 2). Beside the 
Raylcigh-Plateau-like instability we describe a second transversal instability that occurs 
for very small ridges and show that it is related to the spinodal dewetting of a thin film 
on homogeneous substrates. The results of the stability analysis are supported by fully 
nonlinear time simulations. 

The outline of the article is as follows. In Section 11 we introduce the gradient dynamics 
model for the case of a substrate with a wettability pattern. Next, we explain and illustrate 
in Section III C the transversal linear stability analysis procedure employing a sinusoidal 
wettability pattern as example. Our results for a smoothed step-like wettability modulation 
are given in Section 111 B while the cases of convective and diffusive transport are compared 
in Section 111 C. In the subsequent Section IV we discuss numerical time simulations for 
the two different instabilities again considering various different transport mechanisms. We 
conclude and give an outlook in Section V. 



II. THE MODEL 


We employ a mesoscopic continuum model to describe the dynamics of the redistribution 
process under the assumption that the slope of the free surface is everywhere small (long¬ 
wave or small gradient approximation [ODB97, KT07]). In the case of purely convective 
dynamics this assumption allows one to derive an asymptotic model from the governing 
equations of hydrodynamics that describes the time evolution of the height profile of a 
liquid him [ODB97, CM09, ThilO]. It reads 

dkh(x,t) = V- {Q(h)V [P(h, x)]} (1) 

with x = (x, y) T and V = (J^, §~) T ■ Here, Q(h ) = h 3 /( 3//) is the mobility coefficient for a 
fluid where rj is the dynamic viscosity. The generalized pressure P(h, x) is given by 

P(h, x) = — 7 A h — Pt(h, x) , (2) 

where 7 is the liquid-gas surface tension in the Laplace pressure term and 17 (h, x) is the 
Derjaguin or disjoining pressure [dG85, SV09, Isrll], For the latter we use [Pis02] 

n (h, x ) = (! + P^( x )) • (3) 

Here the space-dependent term pp(x) takes into account the chemical stripe pattern on the 
substrate. The dimensionless parameter p corresponds to the strength of the wettability 
contrast, whereas the function p(x) is periodic and describes the geometry of the stripe 
pattern. From now on we use non-dimensional variables h, x, y, t in such a way that 
the non-dimensonal parameters 3 fj, 7 , A, B are all equal to one [nota]. The tilde is then 
dropped in the following. The non-dimensional precursor him height is then uniformly 
h p = 1. Multiplying the term 1 + pp(x) to the disjoining pressure keeps h p constant and only 
modulates the equilibrium contact angle as 60 = y/f (1 + pp(x)). Note that 9 0 is the angle 
in long-wave scaling, i.e., a small physical equilibrium contact angle 6 e = e 6 *o corresponds to 
a long-wave contact angle Oq of 0 ( 1 ). 

Finally, we note that the hydrodynamic thin him equation corresponds to a gradient 

dynamics of the underlying free energy (or interface Hamiltonian [BEI + 09]). The pressure 

is given by the variation P = —- of the free energy F — f dx [HV/rl 2 + f(h , x)l, where 

oh 

the hrst term represents capillarity (interfacial energy of free surface) and the second term 
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represents position-dependent wettability (wetting or adhesion energy or binding potential), 
which is connected to the disjoining pressure by 77 (h, x) = —dhf(h, x) [Mit93, ThilO]. 
Note, that the interpretation of the thin film equation (1) as a gradient dynamics on the 
free energy functional allows one to go beyond the case of convective transport where the 
mobility is Q(h) = h 3 / (3r/) in the case without slip at the substrate and Q(h) = h 2 /{2>rf) in 
the case with strong slip at the substrate [MWW05]. Indeed, by interpreting h as adsorption 
normalised by a constant liquid bulk density, Eq. (1) becomes a general kinetic equation for 
the transport of material as driven by interfacial energies. In consequence, this allows one to 
investigate (i) transport via diffusion of the entire adsorbed film (then one uses the diffusive 
mobility Q(h) oc h as in dynamical density functional theory [AR04, ART10]), and as well 
(ii) transport via diffusion of a surface layer on the deposit as in typical solid-on-solid models 
(then one has a constant mobility Q(h) oc 1 [SS86, MPLS10]). Note, that similar mobilities 
have been introduced in a piece-wise thin-film model for the formation and motion of droplets 
on composite (meltable) substrates [YKP07]. A comparative discussion and interpretation 
of result obtained with the various mobilities is provided below. 

III. LINEAR STABILITY OF RIDGE STATES 

A. Sinusoidal wettability pattern 

To illustrate how a continuation method can be employed in the linear stability analysis, 
we first use a simple harmonic function of one independent variable x as wettability pattern 

0 (x) = sinf^Y (4) 

i.e., the substrate pattern is invariant in ^/-direction. The domain size L is chosen equal 
to the period of the stripe pattern L per . First, we determine steady solutions ho(x) of the 
one-dimensional system. Setting dt h = 0 in Eq. (1) and integrating twice leads to 

d 2 x h 0 (x) + n(h 0 ,x) + C — 0. (5) 

The first integration constant corresponds to a net flux into or out of the integration domain 
and therefore vanishes in the case of horizontal substrates without additional lateral driving 
forces. The second integration constant C is the constant pressure inside the film that 
characterizes mechanical equilibrium. In the following, the constant pressure C is used as a 
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continuation parameter that plays the role of a Lagrange multiplier which ensures a constant 
deposit volume (as the density is assumed to be constant, here this is equivalent to mass 
conservation). 

Such a one-dimensional solution can be extended homogeneously in ^/-direction thereby 
creating a steady solution of the two-dimensional thin him equation (1) that is translationally 
invariant in ^/-direction, i.e., all y derivatives vanish. To determine the stability of such a 
solution, we allow for perturbations that depend in an arbitrary manner on x and are 
harmonic in the y direction. We use the ansatz 

h(x, y, t ) = h 0 (x ) + £ h\(x) exp (/3t + i qy) (6) 

where (3 is the growth rate, q is the transversal wavenumber and £ < 1. Introducing this 
ansatz into equation (1) leads to lowest order in e to the linear eigenvalue equation 

Phi = - Q{h 0 )(dl - q 2 ) 

x [(d 2 x - q 2 )hi + (d h n(h 0 , x))hi] 

- (d x Q(ho))d x [{d 2 x - q 2 )h x + ( d h n(h 0 , x))hi] (7) 

that has to be solved for unknown ft and hi(x). We solve the time-independent equa¬ 
tion (5) and the eigenvalue problem (7) in parallel using pseudo-arclength continuation as 
implemented in the continuation toolbox ALITO-O7p [DKK91, DO09]. The specific set of 
equations treated by ALITO-O7p is given in Appendix A. 

The desired results are obtained through a number of continuation runs: We start with 
the trivial solution h 0 (x) — h — const, and hi = 0, set p, ft = 0 and C = —Il(h,x), whereas 
the parameter q is fixed at some small positive value. We continue this solution as we change 
the wettability contrast p. This yields the solution branches depicted in Fig. 3 a). Figure 
3 b) shows three examples of ridge/trench cross-sections according to the labels in a). 

In the next step, we start at a specific point on the stable part of the solution branch 
for positive p and let the growth rate (3 vary. As h\ is initially zero, and ho is independent 
of (3, and the solution does not change with (3. However, when [3 reaches a value of the 
discrete eigenspectrum, a branching point is detected as a solution branch with hi ^ 0 
bifurcates vertically. In the next run one follows this new branch, effectively only changing 
the amplitude of hi at fixed (3. The norm ||hi|] is used as an additional continuation 
parameter and the run is stopped when ||hi|| = 1. 
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X 


FIG. 3: a) Solution branches obtained from Eqn. (5), (7) using pseudo-arclength continuation. The 
solid and dashed parts of the branches represent stable and unstable solutions, respectively. The 
parameters are h = 3 and L per = 50. b) Top: Solutions of Eq. (5) according to the labels in a). 
Solutions 1 and 2 are unstable. In these cases, more liquid is on the low wettable stripe. Solution 
3 represents a stable drop on the more wettable patch. Bottom: The inhomogeneity function g{x) 
(cf. Eq. (4)). The more (less) wettable area corresponds to g(x) ~ —1 (g(x) ~ 1). 

Figure 4 a) shows the first two eigenmodes of Eq. (7). The symmetric eigenmode (solid 
line) corresponds to a so-called varicose mode and is unstable, whereas the antisymmetric 
zigzag eigenmode (dashed line) is stable [BWR92, TBBB03]. In the unstable g-range the 
ridge is linearly unstable w.r.t. a Rayleigh-Plateau instability (cf. Refs. TBBB03, BKHT11 
for related thin film systems). 

Now one fixes the norm | \h\\ \ = 1 and continues the eigenvalue problem in the parameters 
/3 and q to directly compute the dispersion relation /3(q ) for each eigenmode. The results 
for both, varicose and zigzag mode, are shown in Fig. 4 b). The varicose mode has a 
finite band of unstable wavenumbers 0 < q < q c similar to the case of a homogeneous 
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FIG. 4: a) The two most important eigenmodes of Eq. (7). The solid curve is an unstable varicose 
mode, the dashed curve is a stable zigzag mode, (p = 0.5, h = 3). b) Dispersion relations of the 
two critical modes shown in a). The varicose mode (solid line) has a band of finite wavelengths 
with positive growth rates. Therefore the analyzed ridge is linearly unstable w.r.t. this mode, 
whereas the zigzag mode (dashed) is always stable for the system in question. 

substrate. The growth rate at q — 0 is always zero for the varicose mode: /3(q = 0) = 
0 corresponding to mass conservation. The zigzag mode is always stable for the present 
system. On a homogenous substrate /3(q — 0) = 0 reflects the translational invariance in 
x clirection[TBBB03]. In contrast, here this invariance is broken by the stripe wettability 
pattern implying that (3(q = 0) ^ 0. Note, that the zigzag mode might become unstable 
under driving in x direction, e.g., on an incline [TK03]. 

We employ AUTO-O7p in such a way that the maximum of the dispersion relation corre¬ 
sponds to a saddle-node bifurcation of the curve q(/3). With this trick one is able to follow 
the maximum of the dispersion relation when varying some other parameter through a fold 
continuation. We follow the maximum of f3(q) for the varicose mode varying, e.g., h, until 
f3 ~ 0. A final continuation follows this point in the plane spanned by h and p and directly 
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gives the instability threshold. This allows us to obtain the stability diagram for ridges 
presented in Fig. 5. 



0 0.2 0.4 0.6 0.8 1 1.2 1.4 

p 


FIG. 5: Linear stability diagram in the [p, h] plane. The solid line shows the linear stability 
threshold for a homogeneous ridge solution. For parameters below and above the curve, ridges are 
linearly stable and unstable, respectively. 

For large wettability contrasts p, the stability threshold with respect to h increases rnono- 
tonically with p. This is the behavior we expect from the experimental findings. For small 
p, however, this is no longer true. We will discuss this surprising result at the end of the 
next section. 

B. Smoothed-step stripe-like inhomogeneities 

In the previous section we have introduced our methodology using a simple sinusoidal 
wettability modulation. Now we use a spatial modulation of the disjoining pressure that 
realistically reflects the experimentally employed chemical stripe pattern. We use smoothed 
steps employing 

g(x) = tanh y(x A ~ ( -) ) 

\ ^per J J _ 

x tanh y (V -X A ) ~ (y—X) , (8) 

_ s V x-^per J J _ 

Figure 7 c) shows an example of the function g(x). The parameter L per is the spatial 
period as before, l s L per is the width of the transition region between the more wettable 
stripe (MWS) and the less wettable stripe (LWS). That is, l s defines the sharpness of the 
wettability contrast. The MWS (g(x) ~ —1) starts at XAL per and ends at (1 — XA)L per in 
each period. The smaller x A is, the wider is the MWS. 


14 




We start with a domain size equal to the period of the stripe pattern L per = 50 and set 
xa = 0.3 so that the MWS is thinner than the LWS. The parameter l s is set to 0.03. As in 
the previous section, the starting solution is the trivial one h$ = h = 3. We continue this 
solution as p is varied and obtain the solution branch shown in Fig. 6. Since the MWS is 
thinner than the LWS, the symmetry between positive and negative p, which could be seen 
in Fig. 3 a) for sinusoidal wettability patterns, is broken. Figure 7 shows the solutions for 
p = 0.5 (Fig. 7 a)) and p = —0.5 (Fig. 7 b)) , respectively. For positive p we find a larger 
variety of unstable solutions, e.g., configurations with a small drop on the MWS and a larger 
drop on the LWS. 



FIG. 6: Solution branches obtained by a first continuation in the wettability contrast p. In contrast 
to Fig. 3, we use the energy as a solution measure in order to assess the stability of the solutions. 
The solutions corresponding to the labels 1 — 8 are shown in Fig. 7. The parameters are L per = 50, 
l s = 0.03, h = 3 and xa = 0.3. 

As in the previous section, we directly compute the linear stability diagram on the [p, h] 
plane. The result is shown as the solid line in Fig. 8. In contrast to the example with 
the sinusoidal wettability modulation (Fig. 5), the curve has an essentially non-monotonous 
form. In the region 0.43 < p < 0.87, there are four different stability regions that deserve 
further investigation. To this end, we consider the system for a constant p = 0.5. 

For p = 0.5 we observe three critical values of the mean film thickness h, h\ « 1.1, 
h 2 ~ 1.8 and h 3 ~ 2.4. We continue the stable one-dimensional solution h 0 (x) and its 
critical eigenfunction h\{x) in the mean him thickness h. Figure 9 shows the obtained 
profiles ho(x) and eigenmodes h\(x ), respectively. In Fig. 10 a), we give the maximal growth 
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FIG. 7: Steady solutions corresponding to the labels a) 1 — 5 and b) 6 — 8 of Fig. 6 for p = 0.5 
(p = —0.5), L per = 50, l s = 0.03, h = 3 and xa = 0.3. As p is positive (negative), the more 
wettable area corresponds to the region where g{x) ~ —1 ( g(x ) ~ 1). Here, solutions five (a)) and 
eight (b)) are the only stable solutions, c) The inhomogeneity function g(x), given by Eq. (8). 

rates 

Ana* = max/3(g) = /3(g max ) (9) 

<? 

and the corresponding g max in the unstable regions h\ < h < h 2 and h > h 3 . In the stable 
regions h < h\ and h 2 < h < h 3 , /3 max = g max = 0 . 

For very thin films (h < hi), the one-dimensional solutions are not yet droplets. They 
are better described by piecewise flat films with a higher mean him thickness on the MWS 
(cf. solution for h = 1.05 in Fig. 9 a)). The hat him pieces are sufficiently thin to be 
below the instability threshold for spinodal dewetting on respective homogeneous substrates. 
Therefore, the piecewise hat him solutions are stable, also in two dimensions. 

For hi < h < h 2 , the one-dimensional solutions have a pronounced droplet shape (cf. so¬ 
lution for h = 1.5 in Fig. 9 a)). The corresponding two-dimensional ridge solutions are 
not stable, but this instability is not the Rayleigh-Plateau-like instability that leads to the 
formation of bulges. This can be inferred from the shape of the critical eigenfunction h\(x), 
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p 

FIG. 8: Linear stability diagram in the \p, h] plane for steep wettability gradients between the 
stripes. In particular, l s = 0.03 (solid line) and l s = 0.001 (dotted line). For p = 0.5 there are 
three stability thresholds hi, h -2 and /13. 

which is depicted in Fig. 9 b) for different h. For a Rayleigh-Plateau-like instability, h\{x) 
is supposed to have two symmetric maxima corresponding to the broadening of the ridge 
(cf. Fig. 11). 

In the region h\ < h < h, 2 , however, the function h\{x) has only one maximum. Therefore 
this instability is related to the formation of droplets on the MWS and therefore belongs to 
the growth regimes I to III observed in the experiments, which are sketched in Fig. 2. 

In the region h 2 < h < h 3 the critical eigenfunction still has only one maximum, but the 
growth rates are negative for finite wavenumbers (and zero for q — 0). In this region a ridge 
on the MWS is stable because it is no longer possible for the ridge to form droplets on the 
MWS. This can be seen from the critical eigenfunction at h = 2.1 in Fig. 9 b). It becomes 
broader than the MWS. Since in this region it is not yet energetically favorable for the ridge 
to leave the stripe, the ridge stabilizes. This corresponds to the growth regime V seen in 
the experiments. 

For mean film heights above h 3 it becomes energetically favorable for the liquid to cover 
also part of the LWS and the maximal growth rate becomes positive (cf. Fig. 10 (a)). At a 
critical film thickness of h* « 2.5, slightly greater than h 3 ~ 2.4, the critical eigenfunction 
undergoes a shape transition from unimodal to bimodal (see Fig. 11). The resulting instabil- 
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FIG. 9: Top panels: Stationary solutions ho(x) (a)) and critical eigenfunction h\{x) (b)) for p = 0.5 
and different mean film heights h. At 2.5 < h < 2.6, the transition between a unimodal and a 
bimodal shape occurs. Bottom panels: The inhomogeneity function g(x), given by Eq. (8). 




FIG. 10: Maximal growth rates /3max and the corresponding g max , which are defined in Eq. (9), 
against h for a) p = 0.5 and b) p = 0.3. 

ity is the Rayleigh-Plateau-like instability that leads to the formation of bulges. Therefore, 
the region h > h 3 corresponds to the growth regime VI (cf. Fig. 2). 

For small wettability contrasts (p < 0.43) we observe the same shape transition in the 
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FIG. 11: Comparison of transversal instabilities with a unimodal eigenfunction h\(x) for small film 
heights (top, h = 1.5) and a bimodal h\(x) for larger film heights (bottom, h = 3.0). Shown are the 
steady ridge profiles (solid curves) together with the sum of the ridge profile and the eigenfunction 
multiplied by a small constant 5 (dashed curve). In the first case the contact line stays fixed, in 
the second case it is shifted. 

critical eigenmode, but the maximal growth rate remains positive for h above the curve in 
Fig. 8. This is demonstrated in Fig. 10 (b) for p = 0.3. 

As a next step we take a look at the influence of the sharpness of the wettability transition. 
The dotted line in Fig. 8 shows the stability threshold in the [p,h] plane for l s = 0.001, 
i. e. for a sharper wettability transition. In this case, stable ridges are possible for lower 
wettability contrasts p. On the other hand, for p > 0.57, the onset of the Rayleigh-Plateau 
instability is shifted towards lower h. This is probably due to the fact that the effective 
width of the MWS decreases with decreasing l s . Further decreasing l s does not lead to a 
significant change of the diagram, therefore l s = 0.001 can already be seen as the limiting 
case of a true step function. 

Finally, we remark that the instability threshold in Fig. 8 reaches the value h ~ 1.26 at 
p = 0 as is the case for the sinusoidal wettability modulation in Fig. 5 above. This value of 
h ~ 1.26 corresponds, as it should be, to the analytically obtained threshold of the spinodal 
instability of a flat him, which is given by dll/dh = 0. This reinforces the interpretation that 
the first instability for small p and h is a spinodal instability. As one increases p from zero to 
a small finite values while keeping h constant, the one-dimensional hat him solutions change 
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towards more and more pronounced droplet solutions with increasing maximal film heights. 
This explains why the instability threshold with respect to h decreases with increasing p. 

C. Comparison of diffusive and convective transport 

Equation (5), which determines the equilibrium morphologies, does not depend on the 
mobility term Q(h) as long as it is non-zero, whereas in the eigenvalue equation (7) only 
Q(ho) enters. Hence, the mobility has an effect on the absolute values of 6 when (3 ^ 0, but 
does not influence the position of the instability thresholds — 0 in parameter space. This 
is demonstrated in Fig. 12, where we plot, as in Fig. 10 (a), the maximum of the dispersion 
relation /3 max against the mean film height h for convective and diffusive mobility functions 
Q(h) rs./ h 3 and Q(h) ~ h, respectively (note the semi-log scale in the upper panel). Although 
the absolute values are dramatically smaller for the diffusive mobility Q(h) = h, the stability 
thresholds /3 max = 0 remain at the same position, as can be seen more clearly in the bottom 
panel, where the corresponding q max is shown as a function of h. Therefore, we argue that 
the specific mobility that results from the particular dominant transport process(es) does 
neither influence the linear stability of the ridges nor the observed equilibrium morphologies. 
Only the time scales at which the instabilities grow depend on the mobility. This point will 
be further elucidated in the subsequent section where we consider the time evolution in the 
fully nonlinear regime. 

IV. TIME SIMULATIONS 

In the previous section we have shown that there exist two different types of transver¬ 
sal instabilities that can occur in different ranges of values of the wettability contrast p: 
Inspection of the respective eigenfunctions indicates that for small volume a ridge decays 
into droplets that only cover the MWS, whereas for large volumes bulges evolve that also 
cover part of the LWS. Depending on the value of p the two instabilities can be separated 
by a range of stable ridges. In the present section we test whether this finding of the linear 
analysis holds in the fully nonliner regime by performing direct numerical simulations of 
Eq. (1) in the two parameter regimes. 

The numerical simulations are performed employing the alternating direction implicit 
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FIG. 12: a) Maximal growth rates of the transversal instability for p = 0.5 against h. The dotted 
line corresponds to a convective mobility Q(h ) = h 3 . This curve is the same as in Fig. 10 a), 
but in logarithmic scaling. The solid line corresponds to a diffusive mobility term h. b) The 
corresponding fastest growing wavenumbers. The growth rates are significantly larger for the 
convective instability, especially in the bulge formation regime. The fastest growing wavenumbers 
are almost identical for both mobilities. 

(ADI) method. Thereby the spatial and time derivatives are discretised using second-order 
finite differences. This results in a nonlinear system of equations that is solved at each time 
step using Newton’s method. In each Newton iteration a linear system (Jacobian matrix 
of the nonlinear system) is efficiently solved employing the ADI method ( see Ref. WB03, 
LKF12 for more details). 

The first simulation is performed at low volume, for h = 1.16. This value corresponds to 
the maximum of f3 ma x that coincides to the maximum of q ma x (cf. Fig. 10 a)), in order to 
minimize the necessary computation time and domain size. The latter is then selected as 
50 in x- and 100 in ^/-direction. 

In the course of the time evolution we measure the non-dimensional free energy 

F{t) = f dx dy [^(Vh) 2 + f(h, x)], where f(h,x) now stands for the non-dimensional 

wetting potential. Note that the total time derivative of the free energy is given by 

dxdy. As the mobility function Q(h ) is always positive, the fi- 
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nal expression is negative, i.e., F is a Lyapunov functional. As a consequence, linearly 
stable and unstable steady solutions of the thin him equation correspond to local minima 
or saddle points of the free energy functional F, respectively. The resulting dependency of 
energy on time is shown in Fig. 13. In this energy time series obtained with the convective 



FIG. 13: Time series of the energy during a direct numerical simulation of the thin film equation (1) 
with h = 1.16 for transport by surface diffusion (mobility Q(h ) ~ 1, dash-dotted green curve), by 
bulk diffusion ( Q(h ) ~ h, solid blue line), and by convection ( Q(h ) ~ h 3 , dashed red line). The 
configurations that correspond to the four energy plateaus for Q(h ) = h 3 are plotted in Fig. 14. 

mobility function Q(h ) (dashed red curve) one can identify four plateaus of different length 
where the evolution approaches different steady states of the system. The first three corre¬ 
spond to unstable steady profiles, while the final one corresponds to the absolutely stable 
steady state configuration. All four configurations for Q(h ) = h 3 are depicted in Fig. 14. 

The first configuration, which is assumed very quickly starting from a flat him perturbed 
by small-amplitude noise, is a transversally invariant ridge that is identical to the one ob¬ 
tained above through one-dimensional steady state continuation. In Fig. 14 a) and Fig. 15, 
the solutions resulting from the direct numerical simulation (dotted line) and of the one¬ 
dimensional continuation (solid line) are compared. The fact that this unstable steady state 
is assumed for a considerably long time during the time simulation retrospectively justifies 
the approach of the transversal linear stability analysis of the previous section. After a 
certain time the homogeneous ridge develops transversal modulations that grow exponen¬ 
tially and evolve into the three droplet state that corresponds to the second plateau and is 
depicted in Fig. 14 b). Two of the drops shrink and vanish in subsequent coarsening events 
until only one droplet remains (Fig. 14 c) and d)). The borders between the LWS and MWS 


22 










a) t = 10 3 

50 

35- 


15 

0 
0 

c) t = 6 • 10 3 

50 


35 

15 

0 


50 

y 


50 

y 


I 


J 


100 


I 

I 


100 


b) t = 3-10 3 


1.8 

50 

1.6 

35 

1.4 

X 

1.2 

1 

15 

0.8 

0 


0 

d) t = 10 4 

50 

35- 


15 

0 


50 

y 


50 

y 




J 


2.5 
2 

1.5 
1 

0.5 


100 


I 

I 


100 


FIG. 14: Snapshots of a direct time simulation of the thin him equation (1) with h = 1.16, l s = 0.03, 
xa = 0.3, L per = 50 in the case of purely convective transport (mobility Q(h) = /i 3 ). The borders 
of the MWS is indicated by dashed (black) lines. Figure 13 shows the corresponding time series of 
the energy F(t). The four snapshots at times t = 1000, 3000, 6000, 10000 correspond to the four 
energy plateaus in Fig. 13. 


in Fig. 14 are indicated by dashed black lines. From this one can see that the droplets are 
completely located on the MWS, in agreement with the experimental case (I) and with the 
discussion in the previous section. 

To further characterize the departure of the time evolution from the ridge state of the first 
plateau, Fig. 16 (solid line) shows the time evolution of the logarithm of the difference A F 
between the energy F(t) and the reference energy F 0 of the ridge solution h = /'J 0 (x). One 
observes that for the time interval when the system is close to the unstable ridge solution (cf. 
Fig. 14 a) and b)), the energy A F grows exponentially with a growth rate 2/3 max (dashed 
curve) as expected from the transversal linear stability analysis. Here, /3 max = /3(q max ) is 
the growth rate, calculated for the most unstable mode of the transversally invariant ridge 
solution. Indeed, an expansion of the energy F(t) about the ridge solution ho with the 
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FIG. 15: Comparison between solutions that result from the one-dimensional continuation (solid 

line) and from the direct time simulation (where the dotted line gives the profile of a cross-section 

at y = 0 and t = 1000 of the transversally invariant ridge state assumed for a long intervall in the 

time evolution, cf. Fig. 14 a)). 
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FIG. 16: The solid line gives log \ F(t) — Fq |, where Fo is the free energy of the steady transversally 
invariant ridge state ho and F(t ) is the free energy at time t during the time simulation in Fig. 14. 
The dashed curve shows the exponential growth with a growth rate 2 (3 m that one expects from the 
linear stability analysis w.r.t. transversal modes. 



perturbed solution of the form h = h 0 + e h\(x) exp (fit) reads 

r ij p 

F(t) = F 0 + e / dx— h l (x)e /3t 
. Oh u-u. 


+£ 2 / dx—— hi(x) 2 e 2/3t + h. o. t.. 

J 6h h=ho 

As the linear part in £ of this expansion vanishes at h = ho, the nontrivial term of lowest 


order is the quadratic one, explaining the growth rate of 2 (5. 

Note that the transient behaviour of the energy F{t) changes when employing different 
mobility functions Q(h ). The solid bine and green dash-dotted lines in Fig. 13) give the 
result for linear (bulk diffusion) and constant (surface diffusion) mobilities, respectively. 
Remarkably, one clearly sees that the same energy plateaus appear independent of the 
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chosen mobility, however, their relative duration depends on the particular chosen transport 
behaviour. For instance, with the convective mobility after the ridge configuration, the 
evolution visits the three drop solution that is still visible as a shoulder with the linear 
mobility and is only barely visible (when zooming in) with the constant mobility. The two 
drop solution is clearly visited by all evolutions and the final configuration (corresponding 
to the single drop solution) is the same for all mobilities Q{h ) what is expected as F is a 
Lyapunov functional, independent of the mobility. 
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FIG. 17: Snapshots of a direct time simulation of Eq. (1) at h = 1.9, p = 0.5, l s = 0.015, xa = 0.4 
and L per = 100 and a convective mobility Q{h) = h 3 . Panels (a), (b), (c), and (d) show snapshots 
at times t = 10 4 , t = 2-10 4 , t = 3-10 4 and t = 8-10 5 , respectively. The borders of the more 
wettable stripes (MWS) is indicated by dashed (black) lines. 


Next, we analyse the time evolution in the parameter range where, according to the 
linear stability analysis in the previous section, the second type of transversal instability 
should occur that is related to bulge formation. The wettability contrast is kept at p = 
0.5, but a larger initial him height h = 1.9 is chosen. Figure 17 shows four snapshots 
from this time simulation for a convective mobility Q(h) = h 3 , whereas Fig. 18 gives the 
corresponding dependence of the energy on time for convective (~ h 3 , red dashed line) and 
surface diffusion (~ 1, solid blue line) mobilities Q(h)). In addition, Fig. 19 shows a zoom 
into the energy evolution for the final thousand time steps (in logarithmic scale). There, the 
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FIG. 18: Time series of the free energy during the time simulation of Eq (1) with h = 1.9 for 
transport by surface diffusion (mobility Q(h) ~ 1, solid blue line) and by convection ( Q(h ) ~ h 3 , 
dashed red line). The configurations that correspond to last four energy plateaus for Q(h) = h 3 
are given in Fig. 17. 



FIG. 19: A zoom into the final thousand time steps in the energy plot (18) for the time simulation 
of Eq. (1) with h = 1.9 for transport by surface diffusion (mobility Q(h) ~ 1, solid blue line) and 
by convection (Q(h) ~ h 3 , dashed red line). The letters (a)-(d) refer to snapshots from DNS (cf. 
Fig. 17) that correspond to last four energy plateaus for Q(h) = h 3 (right). 

(a)-(d) letters correspond to the snapshots from the time simulations that correspond to the 
various plateaus. As in the previous case, the transversally invariant ridge is assumed quickly 
(Fig. 17 a)) and forms the starting point for the subsequent morphological transitions. Note, 
however, that here the ridge solution is accompanied by a single small drop, situated on the 
LWS. The existence of the state with this drop results in a shoulder in the energy plot at 
about t = 10 4 (cf. Fig. 19). After a short time, the drop shrinks and disappears and the 
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system stays close to the unstable steady transversally invariant ridge (with energy F r ) for 
some time (see the short energy plateau in Fig. 19 for Q(h) = h 3 close to F r ). Eventually 
the ridge develops modulations in ^/-direction (see panel b) of Fig. 17). These modulations 
grow and form two large bulges (Fig. 17 c) and Fig. 19 at F = F ( i). Finally, one of them 
vanishes in a coarsening event and only one bulge remains (panel d) of Fig. 17 and Fig. 19 
at F = Fb ). From the dashed black line, indicating the borders of wettability regime, one 
can see that in this case the solution is not pinned to the MWS and also covers the LWS. 
The emerged configurations resemble the experimentally observed case VI (see Fig. 2.) 

However, for larger initial him heights h the final bulge configuration itself (shown in 
Fig. 17 d)) can be unstable. Fig. 20 gives in panels a)-c) three snapshots from the time 
evolution for h = 2.8 with convective mobility Q(h) = h 3 , while Fig. 20 d) shows the evolu¬ 
tion of the energy F(t) for convective transport (red dashed line) and for surface diffusion 
(blue solid line). As in the previous cases, the evolution visits similar unstable steady states 
until a state with two bulges is formed. However, this state does this time not only evolve 
via a coarsening event towards a single bulge, but at the same time forms a liquid bridge 
to the next MWS (the system is periodic in ^/-direction). The final liquid bridge state is 
shown in Fig. 20 c). However, the occurance of this final configuration strongly depends on 
the domain size and periodicity of the wettability pattern. A detailed investigation of this 
aspect is not within the scope of the current study and will be pursued elsewhere. 


V. CONCLUSION 

We have employed a mesoscopic continuum model to analyse the redistribution stage 
of deposition experiments where organic molecules are deposited by vapour deposition on 
a solid substrate with stripe-like wettability patterns. In particular, we have investigated 
the stability of transversally invariant ridges located on the more wettable parts of a stripe 
pattern by means of a linear stability analysis of steady one-dimensional profiles that cor¬ 
respond to the transversally invariant ridges on two-dimensional substrates. Employing a 
technique based on continuation procedures has allowed us to perform the linear analysis 
very efficiently. In particular, we have studied the influence of the wettability contrast, the 
mean film thickness (ridge volume), and the geometry of the stripe pattern. 

We have found that there exist two different instability modes that on the one hand result 
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FIG. 20: Snapshots from a direct time simulation of the thin film equation (1) with h = 2.8, 
l s = 0.03, p = 0.5, x a = 0.3, and L per = 50. Panels a), b), c) show snapshots at times t = 10 4 , 
t = 5 • 10 4 and t = 1.5-10 6 , respectively. The borders of the more wettable stripes (MWS) is 
indicated by dashed (black) lines. Panel d) shows time series of the energy F(t) for transport by 
surface diffusion (mobility Q(h) ~ 1, solid blue line) and by convection ( Q(h ) ~ h 3 , dashed red 
line). 


at large ridge volume in the formation of bulges that spill from the more wettable stripes 
onto the less wettable bare substrate and on the one hand result at small ridge volume 
in the formation of small droplets located on the more wettable stripes, respectively. The 
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different modes are reflected in qualitatively different eigenmodes that have their maxima in 
the regions of the strongest wettability gradient and at the stripe centre, respectively. They 
can be identified as a Rayleigh-Plateau instability and a surface instability as in spinodal 
dewetting, respectively. 

In the parameter plane spanned by mean him height h (liquid volume) and wettability 
contrast p one sees that a clear distinction between the two modes appears at some wetta¬ 
bility contrast p c . For p > p c two unstable him thickness ranges exist with a stable range 
of him thicknesses in between. The critical p c decreases with increasing steepness of the 
transition in wettability. In the opposite limit of a sinusoidal wettability pattern, i.e., when 
the transition zone is as wide as the period, no such distinction of two instability modes 
exists. 

Further, we have investigated the role of the analysed transversally invariant ridge states 
in the course of the fully nonlinear time evolution starting from a deposited homogeneous 
him of molecules towards the hnal bulge or drop geometry. This has allowed us to assess 
whether a detailed analysis of ridge stability is meaningful in the context of deposition 
experiments with striped substrates. We have found that even linearly unstable ridge states 
form the hrst ’organising center’ of the time evolution. This is because they (as well as 
all other unstable states) represent saddles in the state space of the system. One may say 
that they ’attract’ time evolutions from the majority of directions in the inhnite-dimensional 
space of system states and later ’expel’ them into one or a few directions and with rates that 
both can be determined by a linear stability analysis of the unstable states. In particular, 
the evolution that starts at a flat him with small perturbations is hrst attracted by the ridge 
state (the corresponding time scale is related to both, the most unstable eigenmode of the hat 
him and the slowest stable eigenmode of the ridge state. After ’visiting’ the transversally 
invariant ridge state the evolution is expelled along the single unstable direction in state 
space. It corresponds to the most unstable eigenmode of the ridge state (the corresponding 
time scale is the inverse of the growth rate of this eigenmode). The importance of the 
ridge state and also of the other unstable states (as, e.g., multidrop states) for the time 
evolution can be clearly seen in the formation of ’energy plateaus’ in all the plots that give 
the dependence of the free energy on time. Each plateau represents one of the unstable 
steady state. This observation is also important for other gradient dynamics systems as it 
indicates how important it is to understand the complete solution structure of a system in 
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order to control its time evolution. It is further to expect that the unstable steady states are 
the ones that are most easily stabilized through imposed controls - be it spatial patterning 
or feedback. 

The conclusion in the previous paragraph is strengthened by our finding that the rele¬ 
vance of the unstable states holds for all the different dominant transport mechanisms that 
we have investigated: the same energy plateaus appear independently of the chosen mobil¬ 
ity, however, we have found that their relative duration depends on the particular chosen 
transport behaviour. This agrees with the fact discussed above that the steady states do 
not depend on mobilities, but that the time constants and eigenmodes do depend on them. 
As particular examples for the influence of the transport mechanism of redistribution on 
the evolution pathway, we have compared a standard long-wave hydrodynamic model in the 
case of no slip at the substrate (purely convective transport, cubic mobility) with diffusive 
transport via diffusion in the him bulk (linear mobility) and at the him surface (constant 
mobility). All the resulting long-wave evolution equations are gradient dynamics models 
on the identical underlying free energy (interface Hamiltonian) just with dif ferent mobility 
functions. As a result, one may now reconsider the him height evolution equation and com¬ 
bine thin him hydrodynamics for relatively thick layers and dynamical density functional 
theory for relatively thin layers by combining the corresponding mobilities[notb]. This will 
allow one to investigate the cross-over between different dominant transport mechanisms 
depending on various parameters. 

Based on our investigations we conclude more specifically that the different experimen¬ 
tally observed equilibrium morphologies do not result from the different transport mecha¬ 
nisms in the liquid or solid state as thought before. Instead, we have found that the observed 
hat and bulged deposits result from different balances of interface energies and wetting en¬ 
ergies. The influence of bulk clastic energies may also play a role that shall be further 
investigated in the future. 

Appendix A: Appendix: Implementation in AUTO-O7p 

For the treatment with the continuation toolbox AUTO-O7p, we have to transform 
Eqs. (5) and (7) into a system of hrst order ODEs on the interval [0,1]. To this end, 
we hrst define the independent variable £ := x/L with L denoting the physical domain size. 
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Next, we introduce the variables 


MO = MW - h , 

(Al) 

fc\ dh ° 

“ 2<0 = dx 

UX x=L£ 

(A2) 
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Here, h denotes the mean him thickness. With the notation Ui(0 = 

d?ij(^)/d^, the system 

of hrst order ODEs reads 


iii =Lu 2 
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iii = - L[n(h + Ml, u 7 ) + C] 

(A9) 

iis =Lu 4 

(A10) 

ii4 =Lu 5 

(All) 

il 5 =Lu 5 

(A12) 

u 6 =L^ +q 2 u 5 d 2 x (hid h II 0 ) 


[(u e q 2 u 4 + d x {hid h n 0 )] 
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with Q o — Q(u\ + h ) and 

d x (hid h IJ 0 ) = 77" (7 + ui, u 7 )u 2 u 3 + T7'(7 + u u m 7 )m 4 

+ n' x (h + ui,u 7 )u 3 , (A15) 

dl(hid h n 0 ) = II"'(h + Ui, u 7 )u 2 2 u 3 

+ n"(h + Ui,u 7 )(dlh 0 )u 3 
+ 277"(7 + Mi, u 7 )u 2 Ui 
+ IJ' (h + Mi, u 7 )u 3 
+ 2n' x (h + Mi, m 7 )m 4 
+ 277" (7 + Mi, u 7 )u 2 u 3 

+ n 'xxQ 1 + «i, wr)«3 • (A16) 

In the last two equations, primes denote derivatives w. r. t. 7 while the index x means a 
derivative w. r. t. x at constant h. The equation for x = u 7 is necessary because AUTO-O7p 
only allows for autonomous ODEs. 
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